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A  REVIEW  OF  ONE-DIMENSIONAL  OCEANIC  MIXED-LAYER  MODELS 

by 

Jon  Moen 


ABSTRACT 


An  overview  of  the  available  one-dimensional  mixed-layer  models  is  presented. 
A  number  of  specific  models,  developed  from  a  common  set  of  basic  equations, 
are  discussed  and  in ter -compared.  While  many  of  the  models  report  satisfac¬ 
tory  agreement  with  observational  data,  they  depend  strongly  on  the  choice 
of  empirical  parameters.  It  is  pointed  out  that  no  particular  mechanism  has 
emerged  to  fully  explain  observations  of  'slab-like'  (vertically  uniform) 
flow  in  the  mixed  layer  and  so,  despite  the  use  of  slab-flow  in  simplifying 
bulk  model  equations,  its  cause  has  not  been  properly  recognized.  That 
some  extra  mechanism  is  required  can  be  seen  In  the  results  of  the  differ¬ 
ential  models,  which  do  not  predict  slab-like  flow  unless  the  wind  sub¬ 
sides.  A  complete  bibliography  is  provided. 
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INTRODUCTION 


Over  the  past  decade  a  number  of  models  have  been  put  forward  in  an  attempt 
to  explain  the  important  time-dependent  features  of  the  oceanic  mixed  layer. 
They  have  met  with  varying  degrees  of  success  when  tested  against  field 
observations  and  laboratory  experiments. 

An  outline  of  the  various  general  approaches  to  one-dimensional,  upper-ocean 
models  has  been  given  by  Niiler  &  Kraus  [l] .  A  comprehensive  treatment  of 
mixed-layer  forming  mechanisms  was  presented  there  and  these  were  unified 
into  a  generalized  bulk  (vertically  integrated)  model.  A  similar,  although 
less  descriptive,  generalized  bulk  model  was  developed  by  Kim  [2].  Another 
very  comprehensive  review  treating  the  historical  development  by  mechanisms 
rather  than  by  a  series  of  over-lapping  models  has  recently  been  presented 
by  Zilitinkevich,  Chalikov,  &  Resnyansky  [3].  Such  approaches  will  not  be 
followed  in  this  review;  instead,  various  mechanisms  and  approximations 
will  be  discussed  as  they  occur  in  relation  to  each  particular  model.  Com¬ 
parison  between  models  will  be  made  easier  by  use  of  a  consistent  notation 
throughout.  Of  the  many  models  only  ten  have  been  singled  out  for  discus¬ 
sion,  with  the  intention  that  this  choice  should  provide  a  clear  overview 
of  the  development  and  present  state  of  mixed-layer  modelling.  This  does 
not  mean  to  imply,  however,  that  models  excluded  from  this  review  are 
unimportant. 

The  models  to  be  reviewed  fall  fundamentally  into  two  categories:  those 
that  specifically  include  the  effect  of  the  earth’s  rotation,  and  hence 
inertial  oscillations,  and  those  that  do  not.  The  latter  begin  with  the 
model  of  Kraus  &  Turner  [4]  and  assume  that  layer-deepening  occurs  primar¬ 
ily  through  downward  turbulent  transport  of  the  turbulent  energy  imparted 
to  the  surface  by  the  wind.  The  turbulent  energy  available  for  mixing  in 
such  models  is  usually  parameterized  as  a  constant  fraction  of  the  downward 
wind-energy  flux  at  10  m.  The  former  category  begin  with  the  model  due  to 
Pollard,  Rhines,  &  Thompson  [S],  who  proposed  that  the  major  effect  of  the 
wind  acting  on  the  sea  surface  was  to  introduce  inertial  oscillations  into 
the  mixed  layer.  The  effect  of  this  would  be  to  increase  the  velocity  shear 
at  the  mixed-layer  base,  thus  reducing  the  Richardson  number  there  and 
causing  the  instabilities  and  turbulent  entrainment  necessary  for  layer 
deepening.  Such  models  assume  a  'slab-like'  motion  of  the  mixed  layer  with 
mean  velocity  shear  only  within  thin  upper  and  lower  boundary  layers. 

Another  model  due  to  Niiler  [6]  combines  both  the  inertial  oscillation  and 
parameterized  wind-stress  approaches.  This  argues  that  since  mixed  layers 
continue  to  deepen  after  the  initial  one-half  pendulum-day  period  (in  which 
all  of  the  deepening  is  predicted  to  occur  by  the  Pollard,  Rhines  &  Thompson 
model  [5])  a  continuous  and  slow  erosion  process  from  surface,  wind -genera ted 
turbulence  must  also  be  present. 

A  later  model  by  Garwood  [7],  although  including  the  effects  of  rotation, 
represents  a  departure  from  the  earlier  models.  Rather  than  assuming  a 
slab-flow  approximation,  a  more  rigorous  treatment  is  carried  out  in  which 
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the  bulk  flow  is  taken  to  be  the  vertically-averaged  flow.  Shear  production 
of  turbulent  kinetic  energy  within  the  mixed  layer  can  then  be  explicitly 
recognized.  The  source  of  entrainment  energy  consists  entirely  of  surface 
production  and  any  shear  production  within  the  entrainment  zone  is  considered 
to  be  negligible.  The  Richardson  number  is  therefore  not  the  critical  para¬ 
meter  determining  entrainment,  and  the  Important  effect  of  rotation  Is  not  in 
reducing  the  Richardson  number  as  in  the  Pollard,  Rhines  &  Thompson  model  [5jt 
but  rather  in  influencing  the  dissipation  time  scale. 

All  of  the  models  to  be  described  employ  the  same  basic  hydrodynamic  and 
thermodynamic  equations,  which,  when  integrated  vertically,  have  simplified 
bulk  forms.  In  models  employing  the  slab-flow  approximation  the  bulk-energy 
equations  become  yet  further  simplified  through  use  of  the  bulk-momentum 
equation  in  parameterizing  energy  production  at  the  lower  entrainment  inter¬ 
face.  In  fact,  the  magnitude  of  the  production  term,  and  hence  the  model 
performance,  depends  strongly  on  the  slab-flow  approximation.  It  will  become 
clear  in  the  course  of  the  review,  however,  that  although  the  observation  of 
slab-like  motion  is  recognized  by  stipulation  in  bulk  models,  possible 
mechanisms  causing  the  slab  have  not  been  proposed  nor  accounted  for  in  any 
way  in  energy-balance  equations. 


4 


SACLANTCEN  SR-47 


1  THE  BASIC  EQUATIONS 

We  now  proceed  by  writing  down  the  mixed-layer  equations  that  will  form  a 
basis  for  the  description  of  all  the  models  contained  in  this  review. 
Derivations,  when  required,  are  presented  in  Appendix  A. 

For  some  of  the  models  we  do  not  need  to  specify  the  mean  velocity  profile 
within  the  mixed  layer.  For  others  we  consider  a  mixed  layer,  depicted 
schematically  in  Fig.  1 ,  to  be  composed  of  three  sublayers:  a  thin  upper 
shear  layer  of  thickness  d  induced  by  the  wind  blowing  over  the  sea  surface, 
a  thick  middle  layer  of  thickness  h-d  containing  no  shear,  and  a  thin  lower 
shear  layer  of  thickness  6  caused  by  the  entrainment  friction  of  the  mixed 
layer  sliding  over  the  stratified  fluid  below.  In  all  the  bulk  models 
temperature  and  salinity  are  taken  to  be  constant  with  depth  throughout  the 
upper  and  middle  layers  and  all  field  variables  are  taken  to  be  horizontally 
homogeneous.  The  stratified  water  underlying  the  mixed  layer  is  taken  to  be 
at  rest  below  some  level  z  =  -H,  where  H  >  h+S.  We  shall  not  be  concerned 
with  the  effects  of  internal  waves  nor,  apart  from  Denman's  model  (Ch.  3), 
of  horizontal  advection. 

The  vertically  integrated  mean  and  turbulent  mechanical  energy  equations  are 
written  as 

3  _  o  _  all 

It  <KEn<>  *  W^o  •  i  »o@'5>  5z  dz  =  0 

-H 
and 

-  o  _  3U 

(Krt)  +  l  3z  dz  +  ^  + 

~H  0 

0 

+  g  /  w'p1  dz  +  e  =  0, 

-H 

respectively,  where 

KE  *  f  p(U  +  u* )  /  2  dz 
m  _H 


KE!  =  °f  p(u 1 2  +  u'2  +  u'2)  /  2  dz  , 

t  _H  1  2  3 

and  e  is  the  total  rate  of  dissipation  of  turbulent  energy  within  the  water 
column.  The  second  term  in  Eq.  1  represents  the  input  of  mechanical  energy 
supplied  by  the  wind  acting  at  the  sea  surface  and  the  third  term  represents 
the  transfer  of  energy  from  the  mean  to  the  turbulent  flow  by  action  of  the 


(Eq.  1) 


(Eq.  2) 
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Reynolds'  stresses  (u'w1 )  on  the  mean  shear.  The  contribution  to  the  latter 
is  assumed,  in  most  models,  to  be  mainly  from  the  upper  and  lower  shear 
boundary  layers.  The  third  term  of  Eq.  2  represents  the  rate  of  working  by 
the  normal  pressure  fluctuations  in  the  wind  field  flowing  over  the  waves 
together  with  a  surface  energy  flux  due  to  breaking  waves.  The  magnitude  of 
this  term  is  largely  unknown  and  its  parameterization  in  the  various  models 
ought,  therefore,  to  be  treated  with  some  caution.  The  fourth  term  of  Eq.  2 
represents  the  total  vertical  mass  flux  (or  negative  buoyancy  flux)  per  unit 
area  within  the  layer  and,  in  the  absence  of  any  solar  penetration,  is  equal 
to  the  rate  of  change  in  gravitational  potential  energy  PE  within  the  layer, 
given  by 


^  o  gal„ 

sf  TO  -  g  /  (i^V)dz  -  =  0,  (Eq.  3) 

-H 

where  I  is  the  surface  value  of  the  penetrating  component  of  solar  radiation, 

I  =  I  eyz,  a  is  the  thermal  expansion  coefficient,  and  c  is  the  specific  heat 

of  sea  water  at  constant  pressure.  Equation  3  was  obtained  by  multiplication 
of  the  local  mass  conservation  equation 

fMjtSV*!')-0  (Eq.  4) 

by  the  vertical  coordinate  before  integration.  Equation  4  represents  the 
combined  conservation  equations  for  heat  and  salt  and  could,  by  applying  the 
transformation 


be  called  a  buoyancy  conservation  equation.  We  prefer  to  work  with  the 
density  variable,  however,  because  then  it  can  more  easily  be  replaced  by 
PraT  for  comparisons  with  models  that  neglect  salinity  effects  on  density. 

Integration  of  Eq.  4  yields  the  bulk  mass  conservation  equation 

h  ^  +  (w V)0  -  <wV).„  t  |  yi-e^b)  =  0.  (Eq.  5) 

Another  expression  for  the  potential  energy  is  obtained  by  considering  the 
observed  geometry  of  a  deepening  mixed  layer,  which  gives 

3  ijr  ,  ~  \  3h  ,  „h 2  Pq 

SI  +  *gh<po  '  P-h-fi)  3T  +  9“?  3F" 

-gf  Ioe'Yh(h+Y_1)  =  0,  (Eq.  6) 

where3fe  is  the  unit  Heaviside  function,  having  the  value  unity  when 
is  positive  and  zero  otherwise.  The  fourth  term  of  Eq.  6  arises  due  to  the 
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density  changes  below  the  mixed  layer  and  is  derived  from  the  term 
-h 

-1  gz  |f  dz 

-00 

together  with  use  of  Eq.  4  .  This  term,  however,  is  normally  neglected  for 
large  values  of  h. 

It  can  be  shown  (see  App.  A.)  that  the  bottom  boundary  condition  for  the 
mass  flux  is 

W)  =X(p0  -  p_h_4)  If  .  (Eq.  7) 

-h 

Using  this  condition,  Eqs.  5  and  6  combine  to  give  another  useful  form 

It  +  (po  -  -  C9j(wV)0  ♦  4  I0] 

-g^  e'Vh(f  +  Y'1)  =  0  (Eq.  8) 

for  the  changing  potential  energy  of  the  layer. 

Finally,  we  write  the  bulk  momentum  equation  for  the  layer  as 
aR 

at  +  -  *  -  '  lo  *  °>  (Eq-  9> 

where  R  is  the  total  mean  momentum  of  the  layer,  given  by 

0  _ 

R  =  p„  /  u  dz. 

-  o  M  - 


If  a  slab-like  motion  is  assumed,  Eq.  9  can  be  written  in  the  form 

^Hb  +  hd  x  Hb)  +  j'b  ^  =  °'  10^ 

Here  the  second  term  represents  the  effect  of  the  earth's  rotation  and 
the  third  the  bottom  boundary  condition 

(Eq-  "> 

-h  • 

which  simply  describes  the  downward  momentum  flux  across  the  lower  boundary 
as  the  layer  deepens. 
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2  KRAUS  &  TURNER  MODEL  (1967) 

We  begin  with  the  theory  due  to  Kraus  &  Turner  [4],  since  it  is  effectively 
the  forerunner  to  a  host  of  similar  theories  to  follow  it.  The  only  impor¬ 
tant  addition  to  their  theory  was  to  come  later  when  Pollard,  Rhines  & 
Thompson  [5]  included  the  effects  of  inertial  oscillations  in  mixed-layer 
deepening,  but  this  will  be  dealt  with  separately  later. 

By  substitution  of  Eq.  3  into  Eq .  2  to  eliminate  the  total  mass-flux  term, 
and  addition  of  the  mean  Kinetic  energy  equation  (Eq.  1)  we  obtain  the 
mechanical  energy  balance 


{It  ^  -  fl  Jo  } +  k  <KEm  +  ^t> 

+  {  Lw * p 1  +  w1^  ]  -  u0«x^)J  +  e  =  0,  (Eq.  12) 

which  was  written  down  but  not  derived  by  Kraus  &  Turner  in  the  form 

w*  +  G*  -  D*  *  0.  (Eq.  13) 


In  Eq.  13  w*  corresponds  to  the  first  term  of  Eq.  12  and  represents  the 
change  in  potential  energy  due  to  mechanical  stirring  or  free  convection. 

We  see  that  it  is  the  rate  of  change  of  total  potential  energy,  minus  that 

portion,  ^  I  ,  due  to  the  penetrating  component  of  solar  radiation.  G* 
yc 

is  the  kinetic  energy  input  from  the  wind,  some  of  which  will  be  dissipated 
by  the  term  D*.  Since,  from  Eq.  13,  D*  =  -e  we  deduce  that 

G*  =  ^KEm  +  +  [WTPT  +  irwt  ]  '  (Eq-  14) 

o 

Which,  if  any,  of  these  terms  Kraus  &  Turner  intended  to  ignore  is  not 
known,  but  taking  into  account  that  G*  is  set  equal  to  a  constant  in  their 
treatment  it  is  reasonable  to  assume  that  those  contributions  to  Eq.  14 
that  could  vary  with  the  inertial  oscillations  were  implicitly  neglected. 
Also,  on  the  time  scales  considered,  the  turbulent  energy  can  be 

presumed  invariant.  We  will  not  attempt  to  separate  out  these  effects  from 
Eq.  14  since  it  would  be  presumptuous  and,  in  any  case,  this  sort  of 
approach  has  been  carried  out  by  others  (e.g.,  Niiler,  [6])  using  slab 
models  and  so  will  be  discussed  later. 

We  note  that  in  Kraus  &  Turner's  application  of  their  model  to  their 
laboratory  results  [9]  no  mean  flow  was  generated  and  so  the  neglect  of 
inertial  oscillations  from  G*,  though  inadvertent,  was  perfectly  in  order. 
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Elimination  of  the  potential  energy  term  between  Eqs.  12  and  6,  and  neglecting 
the  kinetic  energy  term,  yields 

at2  +  3rgh<po  '  p-h-6>  It 

»  G*-  D*  -  &  IQ  -  g|  Ioe”Yh(h+Y_1 ) ,  (Eq.  15) 

which  is  equivalent  to  Kraus  &  Turner's  equation  10  if  the  term  in  e~Y^  is 
ignored . 

The  model  then  consists  of  bulk  equations  describing  two  regimes.  One 
regime  assumes  mixed-layer  deepening,  in  which  case  the  function  '*  =  1 
and  the  energy  balance  (Eq.  15)  is  integrated  in  time  to  give 


'  h2  - 
po 


t  h(t) 

-  g  /  1 

t=t2  h(t=t2) 


_h_6h  dh  =  /  (G*-D*  -g — -)  dt,  (Eq.  16) 


t=t. 


YC 


and  the  mass  conservation  equation  (Eq.  5),  modified  by  the  entrainment 
boundary  condition  (Eq.  7),  is  integrated  to  give 


0  t=t. 


h(t) 


ol 


/  o_h_6dh  =  /  [-(wV)0  -  ]dt  .  (Eq.  17) 

2)  t=t2 


These  two  equations  require  a  foreknowledge  of  the  initial  values  of 
o  (t2)  and  h(t2)  together  with  the  profile  of  p  below  h(t2).  The 
integrands  on  the  right  of  Eqs.  16  and  17  are  specified  functions  of  time. 


In  the  other  regime  the  mixed  layer  is  not  deepening,  in  which  case  =  0 
and  Eqs.  5  and  15  give 


h ,  -  ~vr)_ 

g[-(5fV)0  - 1 10] 

C5J'  =  /  [-(S’trL  -  ^]dt 

°  t=t,  t=t, 


(Eq.  18) 


(Eq.  19) 


Kraus  &  Turner  investigate  a  simple  case  where  the  dissipation  D*  is  set  to 
zero,  the  penetrating  solar  radiation  ignored  by  setting  y=»,  and  a  constant 
rate  of  working  G*  by  the  wind  is  assumed.  The  surface-heating  term 

_  I0 

-(wT)  -  —  is  then  given  a  saw-tooth  distribution  in  time  and  applied 


10 


SACLANTCEN  SR-47 


to  Eqs.  18,  19,  16  and  1 7 ,_i n  that  order.  These  then  reduce  to  a  simple  set 
of  algebraic  equations  in  pQ  and  h.  Note  that  we  have  switched  to  the 

heating  term  here  by  replacing  p‘  by  p  aT'  in  the  density  flux  term  of 
Eqs.  18  and  19. 

The  saw-tooth  heating  function  chosen  is  of  interest  both  because  it  can  be 
compared  directly  with  earlier  laboratory  experiments  by  Turner  &  Kraus  [9] 
and  because  it  simulates  the  daily  and  seasonal  heating  cycles  quite  well. 
Kraus  &  Turner  report,  for  example,  that  a  sinusoidal  heating  function  that 
they  tried  gave  very  similar  results. 

The  model  predicts  the  following  general  features  of  seasonal  variations 
within  the  mixed  layer: 

1)  The  interface  rises  during  the  heating  cycle,  with  a  minimum  depth 
at  the  time  of  maximum  heating. 

2)  The  time  of  maximum  temperature  lags  behind  the  time  of  maximum 
heating  and  minimum  depth  but  precedes  the  time  of  zero  cooling. 

3)  The  temperature  of  the  layer  is  proportional  to  the  square  of  the 
maximum  rate  of  heating  and  increases  linearly  with  the  period  over  which 
heating  occurs.  It  is  inversely  proportional  to  the  rate  of  kinetic  energy 
input  by  the  wind. 

4)  The  depth  of  the  layer  is  directly  proportional  to  the  mechanical 
stirring  rate  and  inversely  proportional  to  the  heating. 

Comparisons  with  observations  at  Ocean  Weather  Station  (OWS)  PAPA  [ioj  in 
the  northeastern  Pacific  is  made  by  first  setting 


» 


where  u*  is  the  water-friction  velocity  and  u*  the  air-friction  velocity 

W  a 

calculated  from  the  drag  coefficient  and  wind  velocities  at  10  m  using  the 
relation 


Choosing  the  average  values  1.7  x  10"  ,  8  m/s,  and  400  cal/cm2/day  for 
Ci o »  U-|q  and  the  downward  surface-heat  transfer  respectively,  gave  the 

results  h_.  =  26  m  and  Tm:iv  =  9°C,  both  lying  within  the  observed  ranges. 

mi n  max 

Kraus  &  Turner  note,  however,  that  their  model  predictions  lie  well  outside 
the  observed  ranges  if  penetrative  radiation  and  hence  convective  mixing  are 
included.  They  argued,  though,  that  dissipation  should  also  be  included  and 
predictions  might  then  be  improved. 
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3  DENMAN  MODEL  (1973) 

An  approach  similar  to  that  of  Kraus  &  Turner  but  using  somewhat  different 
parameterization  and  allowing  for  a  variable  wind  stress  was  adopted  by 
Denman  [8j.  The  purpose  and  design  of  the  model  was  to  obtain  quantitative 
results  over  a  time  scale  of  several  days  rather  than  months  (as  in  the 
Kraus  &  Turner  model),  to  allow  for  a  meteorological  input  varying  in  time 
and  to  include  insolation  and  mean  upwelling  effects  below  the  lower  inter¬ 
face. 

The  mean-flow  kinetic  energy  balance  is  not  in  itself  considered,  but 
instead,  and  in  common  with  other  treatments  to  follow,  the  starting  point 
is  the  turbulent  kinetic  energy  balance.  However,  this  distinction  is  of 
little  significance,  since  the  mean  and  turbulent  energy  equations  are  not 
mutually  exclusive  descriptions  of  the  flow. 

Elimination  of  the. potential  energy  term  between  Eqs.  3  and  6  yields  an 
expression  for  the  total  mass  flux  within  the  layer  given  by  the  third  term 
of  Eq.  20.  The  turbulent  energy  balance  (Eq.  2)  is  then  written  as 

It  <Wt>  +  |  5z  dz  +  ^  +  ^  ]z=0 


j  >  3h  h2  3po  ga  T 

+  {-We0  ’  PW  9T  '  9  T  JT  -  cy  lo 

+  jp.  |+  e  =  0.  (Eq.  20) 

The  second  term  of_Eq.  20  is  labelled  as  -G*,  the  dissipation  e  as  D*. 

3po 

Substituting  for  by  use  of  the  mass  conservation  equation  (Eq.  5)  and 

3(K?t) 

boundary  condition  (Eq.  7),  neglecting  — ^ —  ,  and  rearranging,  gives 


*  (w  +  {£)  =  -1- -  |  2[-G*+D*  -  31  I0(l-e‘Yh)] 

gh(P0-  p  _h_6  ( 

+  gh[(wV)0  +  ~  (l+e“Yh)]|  ,  (Eq.  21) 

where  w  is  a  slow  upwelling  velocity  below  the  well-mixed  layer  included 
to  allow  for  Urge-scale  horizontal  convergence  or  divergence.  Although, 
for  brevity,  w  was  not  included  in  the  basic  equations,  it  is  clear  that 
these  equations  must  describe  deepening  or  shallowing  relative  to  any  mean 
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vertical  motion  of  the  fluid  itself.  The  absolute  derivatives  in  Eqs.  21 
and  22  are  permitted  because  the  only  horizontal  inhomogeneity  is  that 
giving  rise  to  w.  Since  this  is  assumed  to  have  a  very  large  horizontal 
scale  it  can  be  ignored  in  the  equations. 

L 

Next,  by  substitution  for  gh(pQ  -  p_h_5)  ^  *  again  by  use  of  the  mass 
conservation  (Eq.  5),  Eq.  20  becomes 

W-  -  j  G*  -  D*  -  9h(5fV)0  -  gf  loO-Y-1)  -  If  I0e'Yl1  j  •  (Eq.  22) 


The  density  profile  below  the  interface  is  known  at  any  given  time  but, 
unlike  in  the  Kraus  &  Turner  model,  it  is  permitted  to  change  in  time 
according  to  the  insolation  and  upwelling  effects  described  by  the  mass 
conservation 


I  +  “f  *!  vVvz  ■  o. 


(Eq.  23) 


The  density  at  the  descending  interface  can  then  be  calculated  from  the 
relation 


d 

Ht 


0  =  ^ 
p-h-s  at 

dh  3p 

'  M  az 

z=const. 

z=-h-6 

=  -(w 

.  dh.  3p 
+  3t ; '  az 

_  3£SL  j  PYZ 
c  ioe  ’ 

(Eq.  24) 

-h-6 

where  use  of  Eq.  23  has  been  made. 


Equations  21,  22  and  23  are  then  the  ordinary  differential  equations 
representing  Denman's  model.  The  model  allowed  Denman  to  explore  a  variety 
of  wind-  and  heat-dominated  regimes  with  boundary  conditions  varying  in 
time.  An  analytic  solution  is  presented  in  the  simple  case  of  no  heat 
exchanges,  no  dissipation,  no  imposed  upwelling  velocity  w,  a  constant  wind 

speed,  and  a  linear  density  profile  L  =  below  the  mixed  layer.  In  this 

case  it  is  found  that  the  mixed-layer  depth  is  given  by 


h  =  (I|G1)V3  tl/3  .  (Eq.  25) 

Equation  25  agrees  well  with  an  earlier  independent  derivation  by  Kato  & 
Phillips  [ii]  and  also  with  their  laboratory  experiments  and  those  of  Moore 
&  Long  [12], 
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Following  Turner  [13],  Denman  introduced  an  empirically  determined  constant, 
m,  specifying  the  fraction  of  the  downward  wind  energy  flux  at  10  m 
available  for  mixing.  Equating  this  fraction  to  the  mixing  energy  rate 
gives 

G*  -  D*  =  m)T]U10  .  (Eq.  26) 


Next,  using  the  argument  that  the  stress  |t;|  above  and  below  a  thin  molecu¬ 
lar  surface  boundary  layer  must  be  equal,  and  assuming  that  the  mixing 
energy  rate  in  the  water  can  also  be  given  by 

u*  111  '=  ua(pa/po^  111 

=  C10  Ul0(pa/po>*  HI*  (*>•  27> 

where  the  water-friction  velocity  u*  has  been  replaced  by  u*(p  and 

Wi  a  a  O 

the  air-friction  velocity  u*  by  C|q  U^q,  Eqs.  26  and  27  show  that 
m  =  (C10  pa/po^  •  (Eq‘  28) 


ii 


j- 


The  value  of  m  =  0.0012,  estimated  from  Eq.  28  was  used  by  Kraus  &  Turner 
[4]  and  also  by  Denman  [8]  and  Denman  &  Miyake  [14]  in  their  model  compari¬ 
sons  with  observation  at  OWS  PAPA  in  the  northeastern  Pacific.  The  agreement 
with  observation  was  found  to  be  very  good. 

As  mentioned  earlier,  m  should  in  fact  be  an  empirically  determined  constant. 
Neglecting  heat  exchanges  with  the  atmosphere  and  any  radiation  of  energy 
away  from  the  lower  boundary  by  internal  waves,  Eqs.  2  and  3  show  that  the 
rate  of  change  in  the  potential  energy  of  the  layer  is  exactly  equal  to 
that  of  the  available  mixing  energy  -  however  parameterized.  Since  the 
potential  energy  change  in  a  water  column  can  be  calculated  by  taking 
density  profiles  before  and  after  a  storm,  and  since  the  downward  wind  energy 
flux  at  10  m  can  be  integrated  over  the  same  period,  values  of  the  fraction 
m  can  be  estimated.  Note  that  for  this  fraction  to  be  useful  in  a  mixed- 
layer  prediction  model  it  ought  to  be  universal  and  invariable.  We  see, 
however,  that  the  value  of  0.0012  used  by  Denman  &  Miyaka  [14]  does  not 
agree  with  that  of  0.01  estimated  by  Turner  [13]  using  data  collected  in 
the  southwestern  North  Atlantic.  Clearly,  the  fraction  of  energy  transferred 
from  the  wind  into  mixing  energy  is  not  a  constant  and  would  appear,  at  least 
in  the  two  instances  cited  above,  to  be  a  result  of  different  physical  pro¬ 
cesses.  If  we  consider  the  parameter 

o  3u 

G*  =  -[w'p'  +  w‘KE^  ]Q  -  /  pQ  u'w'  dz  (Eq.  29) 

and  replace  the  integral,  using  the  mean  kinetic  energy  equation  (Eq.  1) 
to  give 

G*  =  -CwTPT  +  +  u^.T  -  L 


4 


(Eq.  30) 
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we  see  that,  apart  from  the  sign,  G*  Is  identical  to  the  form  (Eq.  14)  we 
earlier  attributed  to  Kraus  and  Turner  [4] .  However,  Denman  assumes  that 
the  available  mixing  energy  in  the  layer  (i.e.,  that  remaining  after  dissi¬ 
pation)  can  be  equated  to  the  downward  wind-energy  flux  at  10  m.  He  there¬ 
fore  parameterizes  G*-  D*,  rather  than  simply  G*  as  in  Kraus  &  Turner.  It 
is  possible  that  dissipation  could  vary  with  the  inertial  oscillations  but 
it  seems  unlikely  that  this  would  damp  out  completely  the  inertial  variations 
in  G*.  We  must,  therefore,  draw  a  similar  conclusion  as  for  the  Kraus  & 
Turner  parameterization  and  assume  that  since  G*-  D*  was  set  to  a  constant, 
the  inertial  oscillations  contributing  to  G*,  as  represented  by  Eq.  30, 
were  implicitely  ignored  by  Denman. 

Regarding  the  above  mentioned  anomaly  in  m,  Niiler  [ 6 ;,  whose  model  we  shall 
discuss  later,  argued  that  in  Denman  &  Miyake's  OWS  PAPA  data  [ i'4 ]  the 
deepening  observed  was  not  caused  by  inertial  oscillations  because  the 
starting  depth  for  the  mixed  layer  was  already  greater  than  the  depth  for 
which  such  mixing  could  be  responsible.  He  argued  that  this  was  not  the 
case  in  the  data  used  by  Turner  and  so  could  account  for  the  more  rapid 
deepening  and  resulting  higher  estimates  for  m.  On  the  other  hand,  Martin 
and  Thompson  [15],  whose  model  will  be  discussed  in  Ch.  9,  do  not  find  from 
their  experiments  that  Niiler's  wind-wave  source  of  turbulent  energy  is  in 
fact  necessary  to  explain  observed  mixed-layer  deepening. 
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4  POLLARD,  RHINES,  &  THOMPSON  MODEL  (1973) 

The  main  achievement  of  the  model  due  to  Pollard,  Rhines  &  Thompson  [5], 
henceforth  denoted  PRT,  was  to  isolate  a  new  mechanism  for  deepening  the 
mixed  layer.  In  this  model  it  is  proposed  that  the  mixed  layer  moves  ap¬ 
proximately  as  a  slab,  sliding  over  the  stratified  water  beneath  and  cre¬ 
ating  shear  instability  at  the  interface.  The  main  function  of  the  wind  is 
to  create  mean-flow  energy,  up  to  a  quarter  of  which  then  becomes  available 
for  increasing  the  layer  potential  energy  through  turbulent  mixing  at  the 
lower  entrainment  interface.  The  pattern  of  flow  developing  is  one  of 
strong  inertial  motions  initiated  by  a  high  wind,  which  itself  may  have 
only  a  short  duration.  It  is  implicitely  assumed  that  any  shear  production 
of  energy  within  the  top  boundary  layer  [i.e.  that  assumed  to  be  the  main 
mixing  energy  source  in  the  Kraus  &  Turner  (Ch.  2)  and  Denman  (Ch.  3) 
models]  is  completely  dissipated.  Indeed,  in  their  idealized  model,  PRT 
assume  that  this  is  the  only  dissipation  occurring  within  the  layer. 

Most  of  the  deepening  is  predicted  to  occur  within  the  first  half  inertial 
cycle,  during  which  the  slab  energy  has  built  up  to  its  maximum.  Since  the 
layer  cannot  "unmix",  the  potential  energy,  and  hence  the  layer  depth, 
remains  approximately  constant  thereafter. 

In  terms  of  the  mixed-layer  equations  of  Ch.  1,  the  PRT  model  can  be  devel¬ 
oped  in  the  following  way.  The  momentum  equation  (Eq.  10),  which  has  the 
solution 


Mb  = 


fpoh 


(sin(ft),  cos(ft-l)). 


(Eq.  31) 


is  used  together  with  the  mass  conservation  equation  (Eq.  5)  in  its  time- 
integrated  form  (compare  with  Eq.  17) 


hp„  -  /  P 


ah 


0  0  -h-6  ^ 


dt  +  /  [£  I0(l-e‘Yh)  +  (wVU  dt  =  0  (Eq .32) 


and  the  system  is  closed  either  by  stipulating  that  the  overall  Richardson 
number  is  equal  to  unity  for  marginal  instability 


(P0  "  P-h-fi) 

o  =  -gh  _  6  =  1 

po  %  *ub 


(Eq.  33) 


or  by  assuming  the  energy  balance  (Eq.  12)  to  be  adequately  represented  by 
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We  will  comment  on  both  of  these  'equivalent'  closing  assumptions  shortly, 
but  first,  to  demonstrate  their  equivalence,  we  can  substitute  3(FT)/3t 

from  Eq.  8,  estimate  u^t^  from  Eq.  10,  rewrite  Eq.  34  as 


[p0  +  gh(p0'p-h-6^  It 


[gh(wV)Q  -  32*1  io]  =  0,  (Eq.  35) 


and  then,  neglecting  heat  exchanges  with  the  atmosphere  given  by  the  second 
term,  we  obtain  the  relation  of  Eq.  33. 


The  solution  for  the  layer  depth  is  obtained  by  first  writing 

p-h-6  =  £z  If]  .  .  =  h  -  »  (Eq‘  36^ 

z=-h-o  3  z=-h-o 

where  |£  is  the  picnocline  gradient,  assumed  to  be  constant,  and  N  the 

Brunt-Vaisala  frequency.  Then  substituting  into  Eq.  33  for  u^  and  pQ  , 

obtained  from  Eqs.  31  and  32  respectively,  and  once  again  neglecting  heat 
exchanges  with  the  atmosphere,  gives 


h  =  (!-)!  Kkcoift)l1/4  . 

p  L  f2N2 

o 

For  small  time  (1-cosft)  «  and  Eq.  37  reduces  to 


(Eq.  37) 


(Eq.  38) 


which  is  shown  by  PRT  to  be  in  good  agreement  with  the  laboratory  experi¬ 
ments  of  Kato  &  Phillips  ill]. 

To  show  this  agreement  PRT  first  corrected  Kato  &  Phillips'  experimental 
data  to  account  for  such  effects  as  the  depth  of  the  stirring  grid  below 
the  free  surface,  the  momentum  absorbed  by  the  water  throughout  this  depth, 
and  the  time  taken  for  the  grid  to  reach  full  inertia.  The  discrepancy 
rationalized  by  Kato  &  Phillips  to  account  for  the  original  experimental 
data  is  explained  by  pointing  to  errors  in  the  rationalization  used.  For 
example,  Kato  &  Phillips  assumed  that  the  friction  velocity  u*  could  be 
used  to  represent  the  turbulent  velocity  [u'|,  which  may  not  be  the  case, 
and,  more  seriously,  although  they  used  salt  stratification  in  their  exper¬ 
iment  they  cited  Turner's  [16]  results  for  heat  stratification.  Had  they 

1  /4 

used  the  salt  stratification  results,  an  h  «  t  '  dependence  would  have 
been  deduced  that  clearly  would  have  deviated  from  their  data. 

The  undesirable  aspects  of  the  PRT  model,  apart  from  the  assumption  of 
slab-like  flow  while  applying  a  constant  wind  stress,  are  In  the  assump¬ 
tions  that:  1)  all  the  energy  produced  by  wind-induced  shear  is  dissipat¬ 
ed,  2)  no  dissipation  occurs  at  the  lower  interface  production  zone,  and 
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3)  the  condition  for  marginal  stability  at  the  lower  interface  depends  on 
a  bulk  Richardson  number  using  length  scales  of  the  order  of  the  mixed- 
layer  depth. 

The  first  of  these  assumptions  is  considered  by  Niiler  [6]  to  be  incorrect 
and  in  his  treatment  he  allows  only  a  fixed  fraction  of  this  energy  to  be 
dissipated  while  the  remainder  becomes  available  for  mixing.  Niiler  also 
allows  for  some  dissipation  at  the  lower  interface  (assumption  2  above) 
through  introduction  into  the  mean  momentum  balance  of  a  drag  force  for  the 
inertial  oscillations.  As  may  be  deduced  from  Niiler  &  Kraus  [l] ,  however, 
this  term  probably  has  more  to  do  with  radiation  of  energy  away  from  the 
lower  interface  by  internal  waves  than  with  dissipation  of  energy.  Never¬ 
theless,  Wyatt's  criticism  [17]  of  Niiler's  model  does  not  take  this  drag- 
force  term  into  account.  Wyatt,  instead,  develops  two  variations  on  PRT's 
model.  The  first  parameterizes  into  dissipation  a  fixed  fraction  of  the 
combined  upper-wind-induced  and  lower-interface-production  energies  while 
the  second  makes  the  same  parameterization  but  allows  the  parameter  to  vary 
in  time  according  to  PRT's  Richardson-number  stability  constraint.  These 
approaches  will  be  discussed  in  more  detail  in  Ch.  7. 

The  third  of  the  above  assumptions  is  considered  in  some  detail  by  Wyatt 
who,  after  studying  films  of  an  entrainment  interface  taken  by  Kantha  [18  j , 
came  to  the  conclusion  that  the  appropriate  Richardson  number  for  marginal 
instability  depends  not  on  the  mixed-layer  depth,  but  on  a  length  scale 
that  varies  randomly  in  time,  this  length  scale  being  related  to  the  scale 
of  distortions  caused  by  growing  instabilities  at  the  interface.  It  was 
suggested,  although  measurements  were  not  available  in  support,  that  the 
length  scale  should  be  related  to  the  region  of  velocity  variation.  We 
note,  however,  that  using  smaller  values  of  the  critical  Richardson  number 
defined  by  Eq.  33  as  a  way  of  accounting  for  a  smaller  relevant  length 
scale,  Wyatt  did  not  find  much  variation  in  predicted  mixed-layer  depths 
when  using  her  version  of  the  PRT  model. 
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5  NIILER  MODEL  (1975) 

The  model  due  to  Niiler  [6]  is  effectively  a  combination  of  the  models  due 
to  Kraus  &  Truner  [4]  (or  Denman  [8])  and  PRT  [5].  It  predicts  a  strong 
influence  by  inertial  oscillations  on  the  mixed-layer  depth  over  a  time 

scale  N"1  <  t  <  f~^  .  On  both  smaller  and  larger  time  scales  the  shear 
production  of  turbulent  energy  in  the  surface  wind-wave  driven  layer  domi¬ 
nates  . 

Niiler's  treatment  emphasizes  the  separation  of  the  mixed  layer  into  the 
three  distinctive  sublayers  described  earlier  and  shown  in  Fig.  1 .  By 
making  the  assumptions  that  the  upper  and  lower  shear  layers  are  thin  com¬ 
pared  with  the  mixed-layer  depth  and  that  the  Reynolds  shear  stress  in  the 
upper  layer  and  the  mean  velocity  gradient  in  the  lower  layer  are  approxi¬ 
mately  constant,  the  shear  production  term  in  the  turbulent  energy  balance 
(Eq.  2)  can  be  approximated  by  upper  and  lower  contributions  as  (see  App.  A) 

0  _  -  9u  _  _ 

>o  K"')  •  57  dz  =  -<VV  •  I  -  2°  “b‘“b  !£  .  (Eq.  39) 


The  slab-flow  approximation  has  of  course  been  used  here,  so  that  there  is 
no  contribution  to  shear  production  throughout  the  bulk  of  the  mixed  layer. 
Niiler  does  point  out,  however,  that  if  a  velocity  difference  of  only  O.lm/s 
existed  across  the  mixed  layer,  shear  production  could  indeed  become  impor¬ 
tant  there. 

Using  Eqs.  3  and  6  to  eliminate  the  mass-flux  term,  neglecting  the  time 
rate  of  change  of  turbulent  energy,  and  using  Eq.  39  for  the  shear  produc¬ 
tion,  the  turbulent-energy  balance  Eq.  2  can  be  written 


V”b 


ah 


ah 


—v  on  -to-u/’T  T  \  on  _n 

T~  at  -*9h(p0-p_h_6)  yf  -g? 


3pc 

at 


0 


=  "  [wV+w'KE'^-c  •  (Eq.  40) 

The  solar-heating  contribution  in  e"Y^  has  been  neglected  in  Eq.  40.  If 
now  we  consider  only  temperature  effects  on  density,  write 

9  -  I  -  7  P-h-4]  <*•«> 

for  the  heat  content  of  the  water  column,  and  use  Eq.  36  to  represent 
P-h-6  ’  mass  conservation  Eqs.  5  and  7  take  the  form 
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30 


3h  -i 

W  J 


■-£("IP,>  -  I00-e‘Yh) 

o 

=  -  Q0  ,  say  , 


(Eq.  42) 


where  Q0  represents  the  rate  of  heat  exchange  with  the  atmosphere  across 

the  air/sea  interface.  Then  by  use  of  Eqs.  41,  42,  and  36,  the  turbulent 
energy  balance  of  Eq.  40  takes  the  form 


7f  ■  7  “b’ 


=  {(Uo-Ufa)  •  lo  -  CwTP’  +  ^Vo  ■  e}  -  W 


(Eq.  43) 


The  second  term  on  the  right  of  Eq.  43  is  the  rate  at  which  kinetic  energy 
within  the  top  shear  layer  becomes  available  for  mixing  and,  as  in  Denman 
[8],  Niiler  parameterizes  this  quantity  as  a  constant  fraction  mQ  of  the 

downward  wind-energy  flux  U10|x|  at  10  metres.  We  recall  (see  Eqs.  26 

and  29  that  Denman  [8]  did  not  reduce  this  term  to  its  present  form  before 

parameterization  and  so  implicitly  included  the  contribution  ^  “b  Hb  3T 

from  the  lower  boundary,  which  In  Eq.  43  has  been  separated  out  to  be  dealt 
with  in  another  way.  The  fraction  mQ  was  chosen  by  Niiler  to  correspond 

to  the  true  fraction  scaled  relative  to  the  theoretically  derived  fraction 
1  /2 

(C1Qpa/po)  ,  itself  based  on  the  assumption  (see  Eqs.  27  &  28)  that  the 
available  mixing  energy  is  equal  to  u*|x|  .  Thus 

(v\)  *  I  '  ‘  e 

=  m  u10lll 

=  mo(C10pa/po)1/2  U10|x| 

=  m0  u*|x|  #  (Eq.  44) 

This  alternative  scaling  has  no  significance,  however,  other  than  to  make 
uiQ  an  order-one  quantity. 
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To  close  Eqs.  42  and  43,  Niiler  uses  the  momentum  balance  of  Eq.  10, 
modified  to  include  the  loss  of  momentum  radiated  from  the  bottom  of  the 
layer  by  internal  waves,  as 


a 

al 


(hV  +±x  h Mb  =  ? 


cdI^i 


(Eq.  45) 


where  CD  is  a  drag  coefficient. 


Equations  42,  43,  and  45  then  form  the  basic  equations  for  Niiler's  model. 
The  small-time  (t  <  N~^)  solution 


h  %(12  mQ)1/3 


(Eq.  46) 


was  derived  analytically  by  setting  the  surface  heating  Q  and  the  drag 
coefficient  CQ  equal  to  zero. 


Using  constant  values  of  mQ  ,  Qq  ,  and  CD  ,  numerical  solutions  were 
derived  and  shown  to  display  the  following  general  characteristics: 

1)  For  large  times  of  t  »  f'1  the  layer  depth  was  found  to  approach 
the  value 


2) 

3) 


•'ft) 


3/2 


Do  c 


Qoga 


The  surface  production  of  energy  dominated  the  lower  interface  produc¬ 
tion  not  only  for  times  of  t  <  N~^  but  also  for  t  >  f"* 


If  the  initial  layer  depth  was  greater  than  2  T0^3(pQNf)  ^ 
deepening  was  caused  solely  by  surface  processes. 


4)  The  inertial  oscillations,  through  energy  production  at  the  lower  inter¬ 
face,  caused  rapid  deepening  for  times  N'"*  <  t  <  f”"*  but  for  suffi¬ 
ciently  large  times  (t  »  f "^ )  the  same  eventual  depth  was  found  to 
result  without  inertial  oscillations. 

5)  Without  surface  heating  (Q0  =  0)  the  layer  was  found  to  deepen 
without  bound  according  to  the  small-time  solution  of  Eq.  46. 


The  main  weak  point  of  Niiler's  model,  again  apart  from  the  slab-flow 
assumption  common  to  all  such  models,  is  in  its  treatment  of  mechanical 
dissipation  of  energy.  As  mentioned  earlier,  Wyatt  [17]  has  pointed  out 
that  all  of  the  dissipation  in  Niilers*  model  occurs  within  the  upper  shear 
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layer  and  none  within  the  lower  interface  layer.  In  view  of  parameterization 
(Eq.  44),  however,  her  comment  that  a  value  of  mQ  exceeding  unity  would 

imply  negative  dissipation  appears  to  be  unfounded.  Some  allowance  for 
dissipation  of  the  lower  interface  is  tentatively  made  by  Niiler  in  his  drag 
force  parameter  CD  ,  although  this  term  is  more  fundamentally  concerned 

with  radiation  of  energy  away  from  the  production  zone  by  interval  waves. 

The  importance  of  dissipation  was  later  recognized  by  Niiler  in  the  review 
by  Niiler  &  Kraus  [i] .  There  dissipation  was  divided  into  three  separate 
contributions  arising  as  specified  fractions  of  each  of  the  mixing-energy 
sources:  upper  shear-layer  production,  lower  interface  shear  production  due 
to  the  inertially  oscillating  mean  flow,  and  convective  release  due  to  sur¬ 
face  cooling.  Niiler  &  Kraus  also  proposed  that  these  fractions  be  allowed 
to  vary  in  time  such  that  the  predominating  energy  production  mechanism  be 
brought  into  balance  with  dissipation  as  h  -*■  hM  . 


i 
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6  NON- INTEGRAL  MODELS 

6.1  Mellor  &  Durbin  Model  (1975) 

All  bulk  mixed-layer  models  suffer  from  the  assumption  they  must  make 
regarding  an  imposed  mixed- layer  geometry.  A  mixed  layer  is  thus  not  a 
prediction  of  the  models  but  an  a  priori  assumption.  The  model  due  to 
Mellor  &  Durbin  (henceforth  referred  to  as  MD)  [19]  is  of  great  qualitative 
interest,  therefore,  since  not  being  an  integral  model  it  does  not  need  to 
make  this  imposing  assumption.  It  has  been  used  to  add  authenticity,  for 
example,  to  bulk-model  results  by  Wyatt  [17]  and  by  Martin  &  Thompson 
(henceforth  referred  to  as  MT)  [15],  both  with  qualified  success. 

After  making  corrections  to  their  computer  program,  Wyatt  found  a  certain 
qualitative  agreement  between  the  MD  results  and  her  own  bulk-model  results 
but  found  that  the  mixed-layer  depths  predicted  by  the  former  were  consis¬ 
tently  under  estimated.  MT,  on  the  other  hand,  found  that  over  time  scales 
greater  than  about  one  week,  agreement  with  their  three-layer  bulk  model 
was  good.  We  will  comment  further  on  these  important  comparisons  later  but 
first  it  is  appropriate  to  present  a  brief  summary  of  the  MD  model. 

The  approach  adopted  is  similar  to  that  of  Munk  &  Anderson  [20]  in  that 
the  one-dimensional  momentum  and  heat  equations  are  solved  after  parame¬ 
terizing  the  turbulent  velocity  and  heat  fluctuation  in  the  interior  by  use 
of  eddy  viscosity  coefficients. 

The  equations  are 

3u  .  f _  3u  *] 

3T  *  f  x  H  +  [j“r  '  v  SzJ  *  0  (Eq.  47) 

and  _ 

It  +  w  It  +  It  tTV  "  vt  It  ^  =  0  *  (Eq-  48^ 

where  v  and  Vy  are  molecular  diffusion  coefficients  which,  although  small, 

are  included  for  numerical  reasons.  The  turbulence  terms  are  parameterized 
after  Mellor  &  Yamada  [21]  in  the  form 

_  3U 

-vy  -  lq  Sm  Sz 
and  _  - 

-TV  -  *q  SH  !?  •  <*<•  49> 

where  SM  and  $H  are  complicated  functions  of  the  flux  Richardson  number  and 
universal  empirical  parameters  drawn  from  a  bank  of  neutral  turbulent  flow 
data,  q  is  the  rms  turbulent  velocity  calculated  from  the  steady-state 
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turbulent  energy  equation  (cast  in  a  form  involving  another  empirical 
constant),  and  £  is  a  length  scale  estimated  from 


1  +  Kz/£ 

oo 

f°  qz  dz 


/  q  dz 


(Eq.  50) 


with  m,  yet  another  empirical  constant.  In  fact  MD  used  only  ,  which, 
as  proposed  by  Wyatt,  could  be  a  reason  for  very  little  surface°°production 
of  energy  in  their  model. 

Solutions  to  Eqs.  47  and  48  are  obtained  by  integration  in  time.  At  each 
time  step  the  turbulent  fluxes  are  calculated  by  iteration  between  Eq.  49, 

Eq.  50  and  the  expressions  for  S^,  S^,  the  flux  Richardson  number,  and  the 
turbulent  energy  equation  in  q. 

Although  a  large  number  of  empirically  derived  constants  are  required  by 
the  above  approach,  these  constants  are  fairly  well  established  by  previous 
experimental  work.  A  question  remains,  however,  as  to  their  validity  in 
the  open  ocean.  Wyatt,  for  example,  believes  that  certain  of  the  empirical 
parameters  in  the  calculations  of  and  $N  are  not  valid  in  the  ocean. 

Results  from  the  model  show  that  the  mean  velocity  does  not  become  uniform 
with  depth  unless  the  wind  is  turned  off.  In  that  case  the  velocity  shear 
disappears  within  0.2  of  an  inertial  period.  The  layer  was  found  to  deepen 
most  rapidly  during  the  first  inertial  period  and  continue  to  deepen, 
though  progressively  less  rapidly,  for  remaining  oscillations.  Quite  good 
agreement  was  reported  between  the  model  output  and  OWS  PAPA  [14,  22]. 
However,  in  view  of  MD's  claim  not  to  have  adjusted  constants  in  their 
model  to  obtain  better  agreement  with  the  data  and  since  their  model 
results  are  quantitatively  incorrect  by  up  to  30%  because  of  a  programming 
error  pointed  out  by  Wyatt,  doubt  must  be  cast  on  the  initial  choice  of 
constants  and  parameters,  including  possibly  the  surface-drag  coefficient. 
Wyatt  found  best  quantitative  agreement  between  her  model  (to  be  summarized 
in  Ch.  7)  and  that  of  MD  if  a  rather  low  critical  Richardson  number  ( Ri =0 . 25) 
was  chosen  to  represent  marginal  stability.  In  comparisons  with  OASIN  data 
[23]  collected  in  the  Western  N.  Atlantic,  Wyatt  found  that  Richardson 
numbers  did  not  fall  below  1  and  so  the  MD  model  consistently  under¬ 
predicted  the  mixed-layer  depth.  Correction  of  the  earlier  mentioned  length 
scale  from  £  to  £  (Eq.  50),  thus  increasing  the  surface  shear  production 
of  energy,  d?d  not  apparently  improve  the  predictions. 

Martin  and  Thompson's  (to  be  summarized  in  Ch.  8)  comparison  of  their  model 
with  that  of  MD  showed  quite  good  agreement  between  the  predicted  mixed- 
layer  depths,  provided  that  a  sufficiently  long  integration  time  and  strong 
initial  stratification  were  used.  The  two  predictions  were  found  to  con¬ 
verge  to  the  same  value  after  about  one  week.  It  is  worth  pointing  out, 
perhaps,  that  the  comparisons  carried  out  by  Wyatt  were  for  shorter  inte¬ 
gration  times  (up  to  only  three  inertial  periods)  and  did  not  allow  for 
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surface  heat  exchanges.  It  is  also  of  interest  that  MT  did  not  allow  for 
any  surface  production  of  the  mixing  energy  that  Wyatt  had  shown  to  be 
important  in  explaining  JASIN  data. 

The  good  agreement  found  by  MT  between  the  MD  model  and  ensemble-averaged 
data  taken  from  OWS  NOVEMBER  in  the  north  Pacific  was  for  seasonal  varia¬ 
tions  —  again  on  a  much  larger  time  scale  than  that  considered  by  Wyatt. 

Before  returning  to  discuss  a  number  of  other  bulk  models,  mention  must  be 
made  of  two  more  non-integral  models. 


6.2  Marchuk  et  a!  Model  (1977) 

The  approach  taken  by  Marchuk  et  al  [24]  can  be  distinguished  from  that  of 
MD  in  its  more  complete  treatment  of  the  dynamic  turbulence  equations. 

They  do  not,  for  example,  parameterize  the  turbulent  viscosity  coefficients 
in  terms  of  the  Richardson  number  but  rather  solve  for  them  by  using  the 
unsteady  turbulent  energy  generation  and  decay  equations  closed  hy  use  of 
an  empirical  formula.  Marchuk  et  al  tested  their  model  against  the  ob¬ 
servational  data  of  Halpern  [25]  and  found  good  agreement  between  predicted 
and  measured  parameters,  mentioning  in  particular  layer  depth,  current  shear 
across  the  transition  zone,  changes  in  current  due  to  storm,  estimated  tur¬ 
bulent  viscosity  coefficients,  and  estimated  Richardson  numbers.  Examining 
the  velocity  profiles  in  Marchuk  et  al,  however,  shows  a  somewhat  large 
vertical  shear  within  the  mixed  layer,  both  before  and  during  the  passage 
of  a  storm.  This  is  not  consistent  with  bulk-model  assumptions. 


6.3  Warn-Varnas  and  Piascek  Model  (1979) 

The  model  by  Warn-Varnas  &  Piascek  [26]  is  an  example  of  a  higher-order  tur¬ 
bulence  closure  technique  described  by  Mel  lor  &  Yamada  [21] .  Rather  than 
parameterizing  the  effects  of  the  triple  correlation  terms  as  has  been  done 
in  earlier  models  (e.g.  Garwood,  [7])  the  terms  become  variables  in  a 
higher-order  set  of  equations.  In  these  equations  fourth-order  correlation 
terms  are  parameterized  using  various  empirical  constants.  In  fact  this 
approach  involves  the  use  of  a  rather  larger  number  of  empirically  deter¬ 
mined  constants  than  is  needed  in  simpler  models. 

Determination  of  the  relevant  length  scale,  which  is  important  in  modelling 
the  higher-order  correlation  terms,  requires  the  use  of  yet  another  empiri¬ 
cal  constant. 

Running  their  model,  Warn-Varnas  and  Piascek  show  the  mixed-layer  depth  to 
be  three  times  greater  when  including  the  triple  correlation  terms,  yet 
when  comparing  their  model  with  that  of  Mel  lor  &  Durbin  [19]  their  mixed- 
layer  depths  were  found  to  be  only  10%  greater.  They  suggest  that  the 
much  larger  length  scale  chosen  by  Mel  lor  and  Derbin  served  co  compensate 
for  the  deepening  effect  of  the  triple  correlation  terms  neglected. 
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7  WYATT  MODEL  (1976) 

The  model  most  preferred  by  Wyatt  [17]  in  her  extensive  comparisons  between 
existing  models  and  those  modified  and  improved  by  herself  is  the  one  she 
refers  to  as  PW.  We  will  therefore  discuss  only  this  model  here.  The  model 
derivation  will  not,  however,  be  presented  here  since  PW  is  effectively  a 
combination  of  the  models  due  to  Denman  [8]  and  Pollard  et  al  [5],  both 
presented  earlier.  Because  of  this  combination,  the  model  is  closely  re¬ 
lated  to,  and  so  is  to  be  compared  with,  that  of  Niiler  [6],  also  presented 
earlier.  The  model  carefully  distinguishes  between  two  lower-interface 
production  mechanisms:  shear-flow  instability,  modelled  after  PRT  using  a 
bulk  Richardson  number  criterion,  and  steady  shear  production  of  turbulence 
modelled  by  use  of  the  turbulent  kinetic  energy  balance. 

Wyatt  argues  that  Niiler's  (1975)  model  does  not  take  into  account  any 
dissipation  of  turbulent  energy  produced  by  the  mean  flow  shear  at  the 
mixed-layer/thermocline  interface.  Instead,  all  of  this  dissipation  is 
implied  to  occur  within  the  top  wind-wave  shear  layer.  It  was  proposed  by 
Wyatt  that  dissipation  should  occur  as  a  constant  fraction  of  turbulent 
production  within  each  of  the  upper  and  lower  shear  layers.  In  this  case, 
and  assuming  a  slab-flow  approximation  for  Denman's  (1973)  model,  the  tur¬ 
bulent  energy  equation  (Eq.  2)  is  written  in  the  form 


§tFE'=  {  V<VHb)  +  ! 


h}> 


(Eq.  51) 


where  the  first  term  was  obtained  by  substituting  Eq.  3  into  Eq.  2,  neglect¬ 
ing  the  solar  input  I  ,  and  the  second  term  by  substituting  into  Eq.  2  the 
upper  and  lower  shearproduction  contributions  given  by  Eq.  39.  The  third 
term  formulates  the  assumption  that  a  fixed  fraction  of  the  production  is 

lost  to  dissipation.  The  turbulent  tendency  term  yj-  and  the  surface 

were  both  neglected  in  Eq.  51. 


kinetic-energy  flux  term  Lw'p 1 


+  ^o 


Next,  as  first  suggested  by  Turner  [13] ,  it  is  assumed  that  the  surface 
production  rate  _Tq* (Uq-u^)  can  parameterized  as  a  constant  fraction  K  of 

the  downward  wind -energy  flux  at  10  m.  Thus  Eq.  51  is  rewritten  as 


=  [mw  kU1Q  +  IP],  (Eq.  52) 

where  IP  is  a  label  for  the  interface  production  term  $  u^.u^  *n 
this  formulation  the  fraction  of  the  wind  energy  flux  at  10  m  actually 
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available  for  mixing  is  equal  to  rr^k.  This  is  to  be  compared  with  the 

parameter  m  of  Eq.  44  in  our  description  of  Niiler's  model.  The  term  IP  in 
Eq.  52  is  calculated  using  the  slab  momentum  balance  equation  (Eq.  10). 

The  energy  balance  (Eq.  52)  actually  defines  a  model  called  DW  by  Wyatt. 
Model  PW  was  obtained  by  applying  DW  until  the  bulk  Richardson  number 


(p0  ”  P-h.x) 

R.  =  -gh 


po  V“b 


(Eq.  53) 


fell  below  a  critical  value  of  1,  say.  In  this  case  the  energy  equation 
(Eq.  52)  was  replaced  by  the  Richardson  number  (Eq.  53)  (set  to  its  critical 
value)  as  the  closing  assumption.  Deepening  was  thus  assumed  to  be  con¬ 
trolled  by  a  stability  parameter  in  this  way  until  Eq.  52  indicated  that 
the  Richardson  number  was  no  longer  critical.  The  layer  was  found  to 
deepen  much  faster  during  the  period  of  marginal  instability  and,  assuming 
that  the  energy  balance  must  remain  satisfied,  calculations  by  substituting 
Eq.  53  in  Eq.  52  showed  that  the  proportion  of  mixing  energy  dissipated 
during  this  period  was  much  higher.  PW,  then,  was  effectively  the  applica¬ 
tion  of  Eq.  52,  but  with  m  varying  in  time  according  to  the  value  of  a 
stability  parameter. 

A  critical  value  of  Ri  =  1  was  chosen  by  PRT  to  represent  marginal  stabili¬ 
ty.  This  value  was  also  used  by  Wyatt,  but  other  lower  critical  values 
were  tried  in  an  attempt  to  recognize  that  h  in  the  definition  of  Eq.  49 
was  not  really  the  appropriate  length  scale.  The  appropriate  length  scale, 
it  was  proposed,  should  be  related  to  the  thickness  of  the  mean-velocity 
shear  layer,  but  since  this  was  unknown  a  more  rigorous  treatment  was  not 
attempted. 


In  model  comparisons  with  JASIN  data  [23]  Wyatt  found  that  the  critical 
Richardson  number  was  consistently  above  1,  and  so  the  tested  model  PW 
effectively  became  the  model  DW.  By  a  suitable  choice  of  and  K  the 


model  DW  was 
earlier,  the 
(Sect.  6.1). 
was  found  to 
Denman  model 


found  to  give  the  best  fit  to  the  data.  Also,  as  mentioned 
model  was  found  to  give  better  predictions  than  the  MD  model 
For  this  particular  data  set,  surface  production  of  energy 
dominate  over  lower-interface  production  to  the  extent  that  a 
(Ch.  3)  would  also  have  been  a  reasonable  approximation. 
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8  GARWOOD  MODEL  (1977) 

In  the  model  by  Garwood  [7],  entrainment  at  the  lower  interface  and  hence 
layer  deepening  is  controlled  by  a  parameter  expressing  the  ratio  of  buoyant 
damping  to  the  convergence  of  turbulent  kinetic  energy  at  the  interface 


_  JL_ 


p  'w 1 


m4  = 


(Eg.  54) 


This  contrasts  with  the  more  commonly  used  gradient  Richardson  number,  which 
expresses  the  ratio  of  buoyant  damping  to  the  production  of  turbulent  ki¬ 
netic  energy  within  the  entrainment  zone.  Garwood  considers  such  production 
to  be  insignificant  when  compared  with  the  local  convergence  of  turbulent 
energy.  Using  the  dimensionally  derived  turbulent-transport  time  scale 

Tg  =  a1h<w,2>  ,  (Eq.  55) 


with  a1  a  constant  of  proportionality  such  that 


)  = 


* 


(Eq.  56) 


where  the  angle  brackets  '<  >'  represent  vertical  averages.  Such 
averages  are  to  be  distinguished  from  the  earlier  defined  bulk  values 
(  )b,  since  in  the  latter  it  was  assumed  that  the  variable  was  approxi¬ 

mately  uniform  throughout  the  mixed  layer.  Combining  Eqs.  54,  55,  and  56 
yields,  for  the  entrainment  mass  flux, 


(Eq.  57) 


Also  required  for  the  model  are  the  usual  bottom  boundary  conditions  or 
specific  jump  equations  (Eqs.  7  and  11)  describing  the  mass  and  momentum 
fluxes  across  the  entrainment  interface.  These  take  the  forms 

-(wV)_h  =  (<P>  -  P-h-s)  |t  (E<1-  58) 

-(^.h  =  (<u>  -  u_h_6)  .  (Eq.  59) 
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We  note  that_some_error  is  likelyjin  using  these  relations,  since  the 
quantities  <p>  -  P_^_5  and  <u>  -  M-h-6  are  assumed  to  represent  the  jumps 

in  density  and  velocity  respectively  across  the  entrainment  zone.  For  the 
density  this  assumption  is  reasonable,  but  since  velocity  shear  is  allowed 
here  to  exist  throughout  the  mixed  layer  the  velocity  jump  across  the 
entrainment  zone  ought  to  be  considerably  less. 

In  common  with  Pollard,  Rhines  &  Thompson  [5]  and  with  Niiler  [6],  the  model 
requires  the  mass-conservation  and  horizontal -momentum  balances  of  Eqs.  5 
and  10,  which  we  rewrite  here  by  use  of  Eqs.  58  and  59  as 

h|f  <P>  +  (<P>  -  P_h_fi)  §£  +  (wV^  +  ^  =  0  (Eq.  60) 

and 

h^r  <u>  +  (<u>  -  u_h_6)  +  (u'W)0  +  h(fx<]J>)  =  0.  (Eq.  61) 


Also  in  common  with  Niiler  [6]  the  problem  is  now  closed  by  modelling  the 
turbulent  kinetic-energy  equation.  Here,  however,  the  treatment  is  novel 
and  more  rigorous  than  in  previous  models.  The  energy  balance  is  first 
divided  into  horizontal  and  vertical  components,  thus  recognizing  that  the 
turbulent  eddies  are  not  isotropic.  Each  term  of  the  balance  is  then  para¬ 
meterized  by  turbulent-field  modelling  techniques  and,  in  the  case  of  dissi¬ 
pation,  by  a  judicious  choice  of  time  scales.  The  integrated,  component, 
turbulent-energy,  balance  equations  can  be  written  as 


„  <u'  .u‘>  3u  ^  „  u' .u' 

9  (h  +  /°  (TV.  ^  dz  +  f°  ^  (w'^=  )  dz 


Jt  V"  T 


-h-6" 


l  I 


o  p,  3u,  3u9  o  , 

1  +  Tx?  dz  +  l  ,  “l  *'dz 


-h-6  OAl 


-h-6 


3u .  3u . 

— —  — 1  dz  =  0 


+  3  {  .  V  3x .  3x, 

-h-6  j  J 


(Eq.  62) 


and 


L(h<^rr>2)+  f  -S-(wTpT)dz  +  /  L (w'3+  5LE-  )  dz 
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(Eq.  63) 


The  shear-production  term  of  Eq.  62  is  modelled  as 


o  _  3u  9  0  _  3u 

/  (u V )•  dz  «  -  u*  6  U  (0)  +  /  u'w'  -5-  dz  ,  (Eq.  64) 
-h-6  92  -h  92 


where  6  U(o)  is  a  measure  of  the  excess  surface  velocity  due  to  shear. 
It  is  then  assumed  that 


6  U(o)  =  m^  u*  ,  (Eq.  65) 

where  m3  is  a  constant  to  be  determined,  and  that  the  integral  on  the 
right-hand  side  of  Eq.  64  is  given  by  (compare  Eq.  39)  the  expression 


-h-6  3h  .  (Eq.  66) 

2  3t 


As  in  Niiler  [ 6 7,  Garwood  augments  this  shear  production  with  the  surface 
flux  of  turbulent  kinetic  energy  and  pressure  perturbations,  as  represented 
by  the  third  term  in  Eq.  62.  The  lower  contribution  of  this  term  is  asso¬ 
ciated  with  the  radiation  of  energy  away  from  the  mixed  layer  by  internal 
waves  and  is  assumed  to  be  unimportant  in  this  model. 


The  surface  flux  of  kinetic  energy  is  thought  to  be  caused  primarily  by 
breaking  surface  waves  and  is  taken  to  be  proportional  to  u£  for  a  fully 

developed  wave  field.  The  constant  m3  in  Eq.  65  is  thus  adjusted  to 

account  for  breaking  surface  waves.  A  problem  arising  in  Garwood's  treat¬ 
ment,  however,  is  that  it  is  not  known  what  proportion  of  the  surface 
energy  flux  makes  a  contribution  in  the  horizontal  energy  balance  (Eq.  62) 
and  what  proportion  makes  a  contribution  in  the  vertical  energy  balance 
(Eq.  63)  .  Garwood's  method  appears  to  attribute  the  whole  of  this  term  to 
the  former  (see  Garwood's  equations  37  and  38). 

The  pressure  redistribution  term  in  Eq.  62  is  rewritten  as 


o 

/ 

-h-6 


0 

f 

-h-6 


p'  3w' 

P0  9Z 


dz 


=  m2<KEt>  (<KEt>-3<w'2>), 

where  the  first  step  makes  use  of  the  continuity  condition 


(Eq. 
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and 
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the  second  represents  a  form  said  to  be  consistent  with  the  work  of  Rotta 
[27]  and  Lumley  &  Khajeh-Nouri  [28].  The  constant  m2  is  a  parameter 
related  to  the  redistribution  length  scale. 

The  net  rate  of  dissipation  is  considered  to  have  two  contributions:  one 
for  shallow  mixed  layers  where  rotation  is  unimportant,  and  one  for  deeper 
mixed  layers  where  rotation  exerts  a  strong  influence.  The  time  scale 
associated  with  the  former  is  assumed  to  be  proportional  to  the  mixed  layer 
depth  divided  by  the  rms  turbulent  velocity. 

t,  =  h<KEt>'*  (Eq.  68) 

and  the  time  scale  associated  with  the  latter  is  simply 

t2  =  f_1  (Eq.  69) 


The  introduction  of  a  rotational  time  scale  for  dissipation  is  a  new  concept 
in  bulk  mixed-layer  models.  It  is  argued  that  since  rotation  turns  the  mean 
shear  direction  with  depth  there  is  an  inseparable  link  between  rotation 
and  local  shear  production  of  energy.  The  assumption  is  then  made  that  the 
local  rate  of  dissipation  is  related  to  the  local  rate  of  production.  The 
convective  and  rotational  time  scales  are  combined  to  give  a  dissipation 
time  scale  t£  such  that 

=  m^  Tj  +  mgT2  ,  ( Eq .  70) 

where  m1  and  m^  are  constants  of  proportionality. 


Next,  by  dimensional  reasoning,  the  net  rate  of  dissipation  is  taken  to 
equal  h<KEt>  t£_1  ,  which  by  Eqs.  68,  69  and  70  gives 

-1 

3/2  — 


o  3u!  3u!  — 

h  ,  v  9x7  ITT  dz  *  h  <Kle  Te 

h-6  j  J 


/ 


m-|<KEt>  +  m^  hf<KEt>  . 


(Eq.  71) 


Combining  the  contributions  of  Eqs.  64,  67  and  71  ,  the  energy  balance  of 
Eq.  62  can  be  written  as 


a 

al 


(h(<KEt>  -  <w,2>)}  -|m3  u*  +  ^  ~  ^  J 


+  (m2  <KEt>  (<KEt>  "  3<wl2>)} 

2  - 3/2  - 

+  y  (m]<KEt>  +  m5hf<KEt>}  =  0  .  (Eq.  72) 
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For  the  vertical  balance,  Eq.  63,  the  potential  energy  term  is  estimated 
either  by  a  double  integration  of  the  mass  flux,  Eq.  4,  or  by  elimination 
of  the  potential  energy  between  Eqs.  3  and  6.  Using  the  former  method  we 
have 


f-h(iTV). 


♦pS-iy  ‘ol’V)  .  (Eq.  73) 

0 

2  3  w1 d 1 

Ignoring  the  flux  convergence  term  in  ^  (w*  +  -p) »  which  was  not 

accounted  for  by  Garwood,  neglecting  the  rotation  redistribution  term,  and 
using  Eqs.  73,  67,  and  71,  yields  for  the  vertical  energy  balance 


W  (h<w'2>)  *  ♦  <^>0  +  I  ywhy)-'} 

—  i  —  —2 

-  {m2  <KEt>  (<KEt>  -  3<w '  > ) } 

+  m]  <KEt>3/2  +  m5hf<KEt>}  =  0  .  (Eq.  74) 


Equations  57  and  58  combine  to  give  a  third  equation 
ah  <w'2>*  <KE.> 

(<P>-P.h_6)  i -  =  0  .  (Eq.  75) 


Equations  72,  74,  and  75,  together  with  the  bulk  mass  and  momentum  Eqs.  60 
and  61,  then  form  a  closed  set  of  differential  equations  in  h,  <p>,  <u>, 
<wT2>,  and  <KEt>.  The  adjustable  parameters  in  the  model  are:  the  surface- 
heat  flux  term  (w'p1 )  ,  the  penetrating  component  of  solar  radiation  I  , 

the  surface  wind  stress  (u/w' )Q  ,  the  density  P.^  below  the  mixed 
layer,  and  the  velocity  £_h_5  below  the  mixed  layer. 

By  estimating  approximate  values  for  the  empirical  constants  m, ,  m2,  - 

m5  and  choosing  hypothetical  values  for  the  adjustable  parameters,  cGarwood 

explores  the  general  behaviour  of  the  model.  The  important  nog-dimensional 
parameters  turn  out  to  be  a  mixed-layer  stability  parameter  H*  expressing 
the  ratio  of  the  effective  buoyancy  flux  to  the  input  of  mechanical  energy 
by  the  wind,  an  interface  stability  parameter  P*  measuring  the  entrainment 
rate,  and  a  parameter  Z*  related  to  the  mixed-layer  depth. 
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These  are  written 

H*  *  [|iy)0  +  I  I0{1-{hv)'1)J/2  m34 

p*  =  5s  (SV).h/2 1"34 

z.  _  !i  hf 

It  is  found  that  the  rate  of  entrainment  P*  decreases  with  the  mixed- 
layer  depth  Z*  and  that  in  the  retreat  mode  (P*  =  0)  the  maximum  sta- 
tability  H*max  is  not  constant  as  in  other  models  but  is  a  function  of  Z* 
(i.e.  layer  depth,  rotation,  and  wind  stress). 

The  ability  of  the  model  to  predict  a  cyclic  steady  state  on  an  annual  scale 
is  demonstrated  by  using  a  sinusoidal  surface-buyoyancy  flux  and  a  constant 
wind  stress.  This  takes  the  form  of  a  closed  loop  in  the  Z*  -  H*  plane. 
Such  a  cyclic  state  is  possible  because  dissipation  for  deeper  mixed  layers 
is  enhanced  by  planetary  rotation,  which  is  assumed  to  influence  the 
turbulence-dissipation  time  scale. 

Another  distinguishing  feature  of  the  model  is  in  its  prediction  of  a  non¬ 
linear  dependence  of  the  entrainment  rate  P*  on  layer  stability  H*. 

Finally,  it  is  suggested  that  plots  in  the  Z*-H*-P*  space  can  be  used  in 
general  as  a  bases  for  comparison  between  the  various  models. 

Although  Garwood's  model  requires  the  use  of  five  empirical  constants  whose 
values  might  be  both  difficult  to  measure  in  the  open  ocean  and  subject  to 
some  variation  from  one  experiment  to  another,  such  a  model,  with  each 
contribution  to  the  generation  and  dissipation  of  mixing  energy  paramete¬ 
rized  separately,  offers  a  great  versatility  in  application.  On  the  other 
hand,  the  sensitivity  of  the  model  to  the  values  of  the  five  constants  has 
not  been  presented  and  if  the  model  is  sensitive  over  the  expected  range  of 
variation  in  the  constants  it  might  be  said  that  this  sort  of  model  is  un¬ 
necessarily  complicated. 
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9  MARTIN  &  THOMPSON  MODEL  (1979) 

The  most  recent  model  to  be  formulated  is  by  Martin  &  Thompson  [15],  here¬ 
after  referred  to  as  MT.  It  is  a  three-layer  upper-ocean  bulk  model  com¬ 
prising  a  mixed  layer,  a  thermocline  layer,  and  a  deep  layer.  The  reason 
for  including  the  lower  layers  was  to  allow  for  horizontal  advection  effects 
in  a  later  installation  of  the  model  as  part  of  a  three-dimensional  numeri¬ 
cal  ocean  model. 

The  mixed  layer  is  assumed  to  move  as  a  slab  and  its  deepening  rate  to  be 
controlled  by  one  of  two  mechanisms.  The  first  is  stratification-limited 
deepening,  in  which  case  the  model  of  PRT  (Ch.  4)  is  used,  with  deepening 
controlled  by  a  stability  parameter  in  the  form  of  a  bulk  Richardson  number 
(see  Eq.  53  or  33).  The  preference  for  depths  predicted  by  the  PRT  model  is 
based  on  their  good  approximation  to  those  predicted  by  numerical  integra¬ 
tion  of  the  Mellor  &  Yamada  level-2  model  [21]  (henceforth  referred  to  as 
the  MYL2  model).  The  MYL2  model  (presented  by  MD)  was  taken  by  MT,  after 
testing  against  data  and  other  theoretical  considerations,  to  be  the  ideal 
mixed-layer  model.  We  note  however,  that  this  view  is  not  shared  by  Wyatt 
[17],  by  Niiler  &  Krauss  [i],  nor  by  Garwood  [7]. 

The  second  mechanism  is  a  heat-flux-limited  regime  where  deepening  is 
controlled  by  a  dimensionless  parameter 

<J>  =  -  (Eq.  76) 

expressing  the  ratio  of  the  Monin  -  Obukov  and  Ekman-layer  scale  lengths, 
where  Q  is  the  surface-heat  flux  minus  the  vertically  integrated  penetra¬ 
ting  component  of  solar  radiation.  The  appropriate  layer  depth  is  expressed 
by  the  form  [18]. 

u* 

h  =  —  h'(<|>)  , 

f 


where  h'  represents  a  function  in  cp  and  is  determined  from  numerical 
experiments  using  the  MYL2  model . 

The  temperature  (density),  depth,  and  mean  flow  velocity  of  each  of  the 
three  layers  in  the  MT  model  are  determined  by  use  of  one  or  other  of  the 
above  deepening  mechanisms  for  the  mixed  layer,  together  with  application 
to  each  layer  of  the  following  bulk-momentum  and  heat-balance  equations: 

r  li-J 
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where  the  subscript  i  refers  to  the  layer  number,  vm  &  are  molecular- 

viscosity  coefficients  included  mainly  for  numerical  reasons,  and  Qq  is 
the  surface  heat  flux.  The  subscripts  i-J  and  i+J  refer  to  the  upper 
and  lower  boundaries  of  the  i^  layer. 

For  momentum  and  heat-balance  calculations  within  the  thermoclir'e  and  deep 
layers,  only  the  vertically  averaged  velocity  and  temperature  are  of  in¬ 
terest.  The  vertical  temperature  distribution  within  the  thermocline  layer 
must  still  be  retained,  however,  since  it  is  important  in  the  calculation 
of  h  for  the  stratification-limited  regime. 

In  a  comparative  test  between  this  model  and  the  MYL2  model  it  was  found 
that  for  the  stratification-limited  regime,  and  particularly  for  weak  strat¬ 
ification,  the  PRT  model  predicted  deeper  mixed  layers,  with  a  convergence 
of  the  two  after  many  days.  To  some  extent  this  is  consistent  with  the 
findings  of  Wyatt  [17]  although  surface  heating  was  not  considered  in  her 
model.  She  did  find,  though,  that  wind  mixing  through  surface  production 
was  important  -  a  mechanism  not  considered  to  be  necessary  in  the  MT  model. 
Had  the  mechanism  been  employed,  as  it  should  have  been  to  be  consistent 
with  the  MYL2  model,  which  did  employ  it,  the  discrepancies  between  the  MT 
and  MYL2  models  would  perhaps  have  been  greater. 

In  the  heat-flux-limited  regime,  agreement  was  somewhat  better,  with  con¬ 
vergence  of  the  two  models  (MT  and  MYL2)  after  about  five  days.  Agreement 
was  less  good  for  low  values  of  surface  heat  flux. 

The  models  were  compared  with  the  north  Pacific  OWS  November  data. 

Predicted  mixed-layer  depths  were  reported  to  be  consistently  less  than 
those  observed  over  the  one-year  period.  The  following  explanations  were 
offered  for  these  differences: 

1)  Ensemble-averaged  winds  were  used  as  input  to  the  model,  thereby  not 
accounting  properly  for  the  non-linear  response  of  layer  depth  to  wind 
speed,  (i.e.  the  layer  deepens  more  rapidly  for  high  winds  but  does  not 
shallow  again  for  light  ones). 

2)  The  ensemble  for  the  wind  data  was  different  from,  although  overlapping, 
the  ensemble  for  the  oceanographic  data  with  which  the  predictions  were 
compared . 

While  these  explanations  appear  acceptable,  they  remain  to  be  tested  in 
some  way  before  the  MT  model,  based  on  the  PRT  model,  can  be  said  to  be  the 
best  available. 
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CONCLUSION 


Each  of  the  models  presented  has  been  developed  from  the  same  set  of  basic 
equations,  using  a  consistent  notation  throughout. 

Three  distinct  mechanisms  for  the  production  of  mixing  energy  appear  to 
have  emerged: 

1)  Surface  production. 

2)  Shear  production  at  the  thermocline  mixed-layer  interface  by  the  action 
of  the  turbulent  Reynold's  stresses  on  the  mean  current  shear. 

3)  Instability  at  the  thermocline  mixed-layer  interface. 

The  surface-production  mechanism  was  first  considered  by  Kraus  &  Turner  [4] 
who,  considering  it  to  be  the  only  energy  source,  parameterized  it  as  a 
fixed  fraction  of  the  downward  wind-energy  flux  at  10  metres.  The  paramete¬ 
rization  of  this  mechanism  has  remained  essentially  the  same  in  later  models 
due  to  Denman  [8],  Niiler  [6],  Wyatt  117],  Garwood  [7],  Niiler  &  Krauss  [1] 
and  others,  although  varying  degrees  of  emphasis  have  been  placed  on  dif¬ 
ferent  contributing  components.  Niiler  &  Kraus,  for  example,  emphasize  the 
surface-energy  flux  term  [w 1 KE. ]  (see  Eq.  2)  with  surface-shear  production 
_  3u 

f  u/w'  •  Tjy  dz  making  a  small  contribution.  All  earlier  treatments, 

however,  reverse  this  emphasis,  with  surface-shear  production  playing  the 
major  role.  Estimates  of  the  parameterization  factor  m  made  from  obser¬ 
vational  data  vary  by  up  to  an  order  of  magnitude,  thus  casting  some  doubt 
on  its  usefulness  in  a  prediction  model.  One  attempt  to  overcome  this 
problem  was  made  by  Wyatt  [17]  who  allowed  m  to  vary  in  time  according  to 
the  value  of  a  stability  parameter. 

Shear  production  of  energy  at  the  thermocline/mixed- layer  interfaced  re- 

3u 

presented  by  the  integral  across  the  interface  of  the  term  u/w'  •  in 

Eq.  2.  It  was  first  treated  explicitly  by  Niiler  [6]  who,  following 
Pollard,  Rhines  &  Thompson  [5],  used  the  momentum  balance  within  the  mixed 
layer.  Indeed,  this  contribution  can  be  estimated  only  by  use  of  the 
momentum  balance,  since  only  in  this  way  can  the  important  effect  of  the 
earth's  rotation  and  hence  inertial  oscillations  be  accounted  for. 

The  third  mechanism  assumes  that  most  of  the  energy  for  deepening  arises 
from  instability  at  the  thermocline/mixed-layer  interface  and  can  be 
modelled  by  specifying  the  value  of  a  stability  parameter  such  that 
marginal  stability  is  maintained.  Instabilities  will  last  for  only  a  rela¬ 
tively  short  time,  with  a  rapid  return  to  a  marginally  stable  state.  The 
first  to  introduce  this  hypothesis  into  a  model  were  Pollard,  Rhines,  & 
Thompson  [5],  They  employed  bulk-momentum  equations  for  the  mixed  layer 
and  closed  the  system  by  specifying  that  the  value  of  the  bulk  Richardson 
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number  be  unity  during  deepening.  This  model,  with  its  widespread  appeal, 
has  been  used  recently  in  a  three-layer  model  by  Martin  &  Thompson  [15]  to 
describe  deepening  in  a  'stratification-limited  regime'.  The  main  difficul¬ 
ty  in  applying  the  mechanims  is  in  deciding  on  a  suitable  value  for  the  bulk 
Richardson  number  or  in  attempting  to  improve  the  theory  by,  instead,  de¬ 
fining  a  local  gradient  Richardson  number.  Also  the  mechanism  does  not 
account  for  the  direct  flux  of  mixing  energy  through  the  sea  surface. 

In  the  model  by  Garwood  [7],  neither  the  second  nor  third  mechanisms  are 
considered  to  be  important  in  comparison  with  the  convergence  at  the  inter¬ 
face  of  the  turbulent-energy  flux  produced  at  the  surface.  According  to 
Garwood  it  is  the  ratio  of  buoyant  damping  to  this  term  that  controls  the 
entrainment  rate  and  hence  layer  deepening.  Another  inovation  due  to 
Garwood  was  to  allow  for  an  enhanced  dissipation  rate  due  to  planetary  ro¬ 
tation.  This  model  was  thus  able  to  reproduce  the  observed  annual  cyclic 
variation. 

A  fourth  mechanism,  mostly  of  importance  over  seasonal  time  scales,  is  pene¬ 
trative  convection.  This  mechanism  has  not  been  discussed  in  this  report. 

Differential  (non-integral)  models  have  been  considered  through  a  description 
of  the  model  due  to  Mellor  &  Durbin  [19]  .  These  have  the  advantage  of  pre¬ 
dicting  the  formation  of  a  mixed  layer  rather  than  pre-supposing  its  exis¬ 
tence  and,  in  particular,  they  do  not  need  to  assume  a  vertically  uniform 
velocity  profile.  Indeed  the  models  due  to  Mellor  &  Durbin  [19],  Marchek 
et  al  [24]  and  Warn-Varnas  &  Piacsek  [26]  ail  indicate  the  development  and 
persistence  of  vertical  shear  throughout  the  entire  mixed  layer,  provided 
that  the  wind  does  not  abate  appreciably.  Since  the  principal  justification 
for  the  integral  or  bulk  model  (apparently  supported  by  observations)  due  to 
Gonella  [29],  Pollard  &  Millard  [30] ,  Pollard  [31]  and  Pollard  et  al  [5]  is 
that  the  velocity  has  approximately  no  vertical  shear  through  a  large  portion 
of  the  mixed  layer,  some  mechanism  to  diminish  the  shear  would  appear  to 
be  missing  from  the  differential  models.  This  same  mechanism  way  well  not 
be  accounted  for  properly  in  the  integral  models  either. 

This  report  does  not  make  extensive  comparison  between  models  and  observa¬ 
tions  but  rather  is  intended  to  clarify  theoretical  differences  between 
the  various  models  themselves. 
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APPENDIX  A 


The  notations  used,  when  not  defined,  are  those  of  the  Main  Text.  Tensor 
notation  is  used  in  some  of  the  equations. 

A.l  The  mixed-layer  momentum  equation 

The  momentum  balance  can  be  expressed  in  the  form 


3u-  3u.j 

"St  +  uj  IxT 
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axi  <p0> 


+  eijk  “j 
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3x -3x . 

J  J 


0  , 


(Eq.  A.l) 


where  is  the  tensor  operator  (e^ jj<  8j  \  *  5.  x  u)  ,  ^  is  the 

earth's  rotation  vector,  p  is  the  pressure  variable,  and  v  the  molecu 
lar-diffusion  coefficient.  Setting  ui  =  u^+u.|  and  taking  the  average  in 

time  gives 
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(Eq.  A. 2) 


Applying  the  horizontal -homogeneity  condition 


ition  =  S*^  =  °) 


and 


neglecting  the  molecular-diffusion  term  leaves,  for  the  horizontal  momentum 
balance. 
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(Eq.  A. 3) 
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Next  we  set  the 
condition 


Sxj 


upwelling  velocity  w  to  zero  and  apply  the  continuity 
0  (Eq.  A. 4) 


to  Eq.  A. 3  to  give 

5z("^)  *  0 

and,  after  vertical  integration  from  below 


7rr(  /  u  dz)  +  f  x  (  /  udz)  +  (w'u1 ) 

O  t  u  ii  O 


-H 


-H 


the  mixed  layer, 
-  (w'u’)_H  =  0. 


(Eq.  A. 5) 


(Eq.  A. 6) 


The  last  term  on  the  left-hand  side  of  Eq.  A. 6  vanishes  if  no  motion  is 
assumed  at  Z  =  -H  .  Equation  A. 6  then  leads  directly  to  the  mixed-layer 
momentum  equation  (Eq.  9)  of  the  Main  Text. 


A. 2  The  lower  momentum  flux-boundary  condition 

Term-by-term  integration  of  Eq.  A. 5  across  the  lower  interface,  remembering 
that  u  =  w  =  0  for  Z  <  -b-6  ,  gives 


-h  3u  a 

'  3Tdz  ’  *?<«<<  ♦«<«’))♦  u  „  &  • 

-h-o 

-h 

/  f  x  u  dz  =  6(f  xu  h)  +  o(62)  , 

-h-6  ~  _  -  -  -n 


{h.{  I?  =  (“^-h  • 

Letting  6  -►  o  leaves 


3h 

‘i-hsT 


(Eq.  A. 7) 


Equation  11  of  the  Main  Text  is  then  obtained  by  assuming  a  slab-flow 
mixed  layer. 
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A. 3  The  mean  flow  kinetic  energy  equation 

Multiplying  Eq.  A. 2  by  u ^ ,  using  the  continuity  condition  of  Eq.  A. 4,  and 
neglecting  molecular  diffusion,  leads  to 


u? 


s  a  _  _  3U,  a  2 

StM  +  ^3<7(ui  UjUT)  “  ujui  9x7  +  uj  3x7  (r" 

J  J  J 


) 


+  _SL  w  p  =0  . 

po 

Applying  horizontal  homogeneity  and  neglecting  upwelling  leaves 
~  il  »u  a  _  _  _ 

|t( - ) +  )±y)  -  yX*az  =  0  •  (Eq.  a.8) 


Vertical  integration  of  Eq.  A.8  from  Z  =  -  H  then  yields  Eq.  1  of  the 
Main  Text. 


A. 4  The  turbulent  kinetic  energy  balance 


The  turbulent  momentum  equation  is  obtained  by  subtracting  Eq.  A. 2  from 
Eq.  A.l  to  give 


3u.'  3u!  3uT  3u.  _  3u! 

W  +  uj  “5*7  ‘  uj  “5x7  +  uj  5x7  +  uj  5x7 

J  J  J  J 


+  s._  sa!  +  a  (El)  .  v 

i3  p_  3x- 'p_ ' 


i  po 


32  ui 

3x -3x .  =  0 

J  J 


Multiplication  of  Eq.  A. 9  by  u^  and  averaging  yields 


5t 


(-1  )  + 


- TTT2 

Ui  - IP"  ,  a  _ 

u  *  "  w  6  > .  u  J  —  +  u  •  u  •  a  u  • 

ij  i  pQ  i  j  3x7  ui 

32 


+  (ui  p,)  +  °j  aly  (“r  }  “  v  a'x^'x"  ("H 


au'j  Su'j 

+  ''3xrLsj r  *° 

J  J 


(Eq.  A. 9) 


(Eq.  A. 10) 
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Applying  horizontal  homogeneity  and  neglecting  upwelling  leaves 


+  r3-  (u|p')  -  v 


82  (u-±) 

3¥s7v2/ 


+  e  =  0 


(Eq.  A. 11) 


which,  after  vertical  integration  from  Z  =  -H  and  neglecting  the  upper  and 
lower  contributions  from  molecular  diffusion  leads  to  Eq.  2  of  the  Main  Text 


A. 5  The  Reynolds  shear  stress  contribution 

Referring  to  Fig.  A.l,  the  integrated  Reynolds  shear  stress  term  can  be 
written  as 

o  9u  o  3u  -h  3u 

/  (u  V )  dz  =  /  (uV)  dz  +  /  (uV)  dz  (Eq.  A. 12) 
-H  -d  -h-6  c 


MIND  - ► 


For  the  first  term  on  the  right-hand  side  of  Eq.  A. 12  a  constant-stress 
layer  is  usually  assumed,  so  that 

3u 


o 


(Eq.  A. 13) 
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For  the  second  term  on  the  right-hand  side  of  Eq.  A.  12  we  let  d  -*■  o  and 
approximate  the  stress  in  the'interval  [-h,  -h-d]  by  boundary  condition 
Eq.  A. 7,  so  that 


-h 

/ 

-h-d 


3u 

(uV)  dz 


-  *  -Is  5  *  S 


dz 


“  ’  H.  3h 
~~2 —  91 


(Eq.  A. 14) 


Expressions  A. 13  and  A. 14  combine  to  form  the  energy  production  term  of 
Eq.  39  of  the  Main  Text. 


A. 6  Potential  energy  for  a  deepening  mixed  layer 

The  changing  potential  energy  of  a  deepening  mixed  layer  can  be  expressed 
by 


-jp  «  Lim  A  (  °  p(t+At)z  dz  -  /  p(t)z  dz)  .  (Eq.  A. 15) 

At  -*0  -a> 

Next  we  consider  the  geometry  of  Fig.  A. 2. 


O 

a 
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and  express  the  difference  between  the  braced  integrals  of  Eq.  A. 15  in  the 
interval  [-h,  o]  as 

f  Ap  z  dz  =  Ap  %  (Eq.  A. 16) 

-h  0  0  c 

and  in  the  interval  [-h-Ah,  -h]  as 

-h 

/  ((p  +  Ap  )  -  Tz)  z  dz  =  -  p  hAh  +  rh2Ah  +  o(Ah2,AhAp  1.  (Eq.  A. 17) 

-h-Ah  0  0  0  o 

Combining  Eqs.  A. 16  and  A. 17  and  substituting  them  into  Eq.  A. 15  yields 

“  =  ‘  9  ?  ~  9h(P0  ‘  rh)  It  +  /  9Z  dz  •  (Eq.  A.18) 

9t 

We  then  set  rh  =  and  note  that  Eq.  A.18  is  equivalent  to  Eq.  6  of 

the  Main  Text. 


A. 7  The  lower  mass  flux-boundary  condition 

Term-by-term  integration  of  the  mass  conservation  equation  of  Eq.  4  of  the 
Main  Text  across  the  lower  interface  gives 


l  .  ll  dz  =  It  (6p  +  +  p-h  It  '  p-h-6  It 


-h 

f 

-h-<5 


I 


o 


letting  6  -*•  o  leaves 


3h 


(w'P').h  -  -<P.h  -  p.h.a)  It 


which  is  the  boundary  condition  in  Eq.  7  of  the  Main  Text. 
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